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Abstract. We investigate the eigenvalue statistics of random Bernoulli 
matrices, where the matrix elements are chosen independently from a binary set 
with equal probability. This is achieved by initiating a discrete random walk 
process over the space of matrices and analysing the induced random motion of 
the eigenvalues - an approach which is similar to Dyson’s Brownian motion model 
but with important modifications. In particular, we show our process is described 
by a Fokker-Planck equation, up to an error margin which vanishes in the limit of 
large matrix dimension. The stationary solution of which corresponds to the joint 
probability density function of certain well-known fixed trace Gaussian ensembles. 


1. Introduction 


Random matrices were first introduced by Wishart in the 1930s in the context of 
multivariate statistical analysis [I]. However it was not until the late 1950s and 
early 1960s that random matrix theory (RMT) was significantly developed by Wigner, 
Dyson, Mehta and others J2[[4]. This interest was initiated by Wigner’s observation 
that eigenvalues of large random matrices and scattering resonances of large atomic 
nuclei shared the same statistical properties j2j . RMT has progressed in a similar 
vain ever since, as the development of new methods has paralleled the application 
to many new areas of mathematics and physics, including, for example, quantum 
chaos, disordered systems, number theory, field theory and financial mathematics [5]. 
One of the main reasons for this outstanding success has been the observation of 
universal features displayed across a wide range of systems, allowing RMT to predict 
the properties of many emergent spectral phenomena, regardless of the intrinsic 
microscopic details of the system. Understanding the reasons for this universality has 
been, and still is, one of the prime goals in RMT and related fields. Much progress has 
been made both in the theoretical physics community (see e.g. [5j|7 and references 
therein) and in mathematics (see e.g. |5][8] and references therein). 

Recently, significant progress has also been made in proving universality for 
Wigner matrices - commonly referred to as the Wigner-Dyson-Mehta conjecture (p.7 
of [9], see also 10,11 ). In essence, this states that for large matrices, in which the 
elements are independent random variables, the distribution of eigenvalues should 
converge to the well known Gaussian ensembles, regardless of the precise distribution 
of the matrix elements. Two complementary approaches have been developed in this 
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respect. The first, based upon the initial ideas of Johansson 12 , is the so-called 


heat-flow method, which analyses the stochastic motion of eigenvalues under a small 
Gaussian convolution 13-16 . The other compares the eigenvalue statistics under the 


element-wise swapping of distributions in a random matrix 17 19 


Remarkably, this universality extends to Wigner matrices in which the entries 
are randomly chosen from a set consisting of just two elements 17 20 , so-called 
Bernoulli ensembles. The methods have also been adapted to treat sparse random 
matrices and distributions with non-zero mean 21 22 . However all these results 


rely on the independence of the matrix elements and a method for relaxing this 
condition for Bernoulli matrices is thus far out of reach. Although universality results 
have been shown in other classes of correlated matrices using orthogonal polynomial 


techniques 23-25 


One of the principle reasons for studying the spectral properties of Bernoulli matrices 
is their inherent relation to graphs, or networks, since this binary choice can be 
used to represent the existence/absence of a connection between vertices, or nodes. 
These random graphs are used extensively in areas such as computer science, network 
modelling, statistical physics and neural networks to name but a few and analysis is 
often based upon their spectral properties. Therefore, a comprehensive understanding 
of this universality is of central importance. Often this binary choice is given by the 
set {0,1}, in which case the matrix is referred to as the adjacency matrix of the graph. 

By imposing different constraints on the matrices one can create ensembles of 
random graphs, that, by extension, can serve to characterise these Bernoulli ensembles. 
Some examples of adjacency matrices associated to well-studied types of graphs are 
listed below. 


(i) Undirected graphs are given by symmetric adjacency matrices A in which the 
elements A, 7 = Aji are equal to 1 if vertices i and j are connected and 0 otherwise. 

(ii) Directed graphs occur when the connection from vertices i to j is independent of 
the connection from j to i. In this case Aij is 1 if there is a directed edge from i 
to j and 0 otherwise. 

(iii) Tournament graphs are unidirectional graphs, meaning all pairs of vertices (i,j) 
are connected with an edge admitting a certain orientation. The matrix elements 
thus satisfy A 3i = 1 — A^ 3 £ {0,1}. 

(iv) d-regular graphs are undirected graphs in which every vertex has exactly the same 
degree, or valency, equal to d. Therefore each A satisfies / A i3 = d, V j. 

(v) Regular tournament graphs are tournament graphs with an odd number of vertices 
N, in which each vertex i has an equal number, N ~ l , of directed edges pointing 
towards and away from i. In other words yv A 33 = jv ~ 1 , V i. 

An ensemble of random graphs is then created by endowing a probability measure 
on the set of all graphs. The most common choice is the uniform distribution, so that 
each graph is equally likely to be selected. In the first three examples above, this 
equates to having iid random variables for the matrix elements in which 0 and 1 are 
chosen with probability 1/2, without any further constraint. We shall therefore refer 
to these ensembles as unconstrained and those with non-independent matrix elements, 
such as random d-regular graphs and random regular tournaments, as constrained. In 
the present manuscript we shall be concerned solely with the former. 
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Returning to the adjacency matrices, then all the ensembles outlined above have 
been shown, 21 22 , or are conjectured 26 27 , to display local eigenvalue statistics 


consistent with RMT, after appropriate rescaling. However these matrices possess two 
properties which render this comparison less than straightforward. The first is that 
E[Tjj] ^ 0, in contrast with the Gaussian case in which the mean is zero. The second 
is that by the Perron-Frobenious theorem, the largest eigenvalue of the adjacency 
matrices is often clearly separated from the rest of the spectrum. We therefore define 
for all A 


B = 2A~ J + I (1) 

where J denotes the matrix in which every element is 1 and / is the identity. B is the 
“balanced” counterpart of A and it does not suffer from the deficiencies mentioned 
above. But for a trivial shift and a factor 2, the difference between A and B is a 
matrix of rank 1, which ensures that the spectral properties of the A and B ensembles 
are closely related. 


In a series of papers, of which this is the first, we introduce a random walk approach 
to the spectral statistics of Bernoulli ensembles. The main purpose of the present 
paper is to introduce the approach and apply it to various unconstrained ensembles. 
We are inspired by Dyson’s Brownian motion model 28 , which returns the joint 


probability density function (JPDF) for the eigenvalues of the canonical Gaussian 
ensembles as the stationary solution of an induced Fokker-Planck equation. However 
Dyson’s model has a particular drawback - it is inextricably linked to the Gaussianity 
of the matrix elements. Here we overcome this restriction, obtaining the same Fokker- 
Planck equation as Dyson, but at the expense of an error caused by having a discrete 
time step dependent on the matrix size N. Therefore, whilst Dyson’s model works 
exactly for any matrix size N > 2, our approach only becomes valid when N 1. 
This relies on taking a series of discrete time steps and showing the linearity in time 
of the first two moments of the evolution kernel - as one would expect from arguments 
pertaining to the central limit theorem. 

The paper is organised as follows: In Section [2] we set up the discrete random 
walk process in the space of matrices and show it is equivalent to a random walk on 
the hypercube. The latter has been well studied previously 29 ■ 31 and is known to be 
ergodic with a uniform stationary distribution. We are interested in how this process 
induces a random walk in a restricted space of variables (in our case the spectrum) 
and the transition from a discrete to a limiting continuous space-time description. 
Therefore, for the purposes of clarity, we first provide a simple example of this limiting 
transition in Section [3j This is the time dependent probability distribution of the 
walker to be at a given distance from the origin 29-31 , which converges to the 


familiar (continuous) Ornstein-Uhlenbeck process in the limit of large hypercubes. 
This shows how entropic forces may arise from a simple reduction of variables and the 
same ideology will be used in order to investigate the spectral statistics of Bernoulli 
matrix ensembles in the remainder of the paper. 

Our main result is outlined in Section [4j Here, using the method of Kolmogorov 
(see also [33) ), we study the evolution of observables under the random walk 
in the space of spectra. This leads to a Fokker-Planck equation with an error that 
decreases like 0(N~ e ), with a value e = 1 — (2a + 7 + c) > 0 that comes from 
choosing our observable to have a spectral subset of size 0(N 7 ) and a resolution 
0(N~A+ a )y The final contribution comes from viewing the random walk in a time 
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window At = 0(N C ). Thus, for large N the use of the continuous Fokker-Planck 
equation suffices to describe the spectral evolution with a resolution better than the 
mean level spacing. In subsequent sections we obtain the precise form of this Fokker- 
Planck equation for three different examples of Bernoulli ensembles, achieved through 
a careful analysis of the moments of the transition probability. The resulting equations 
are exactly those obtained by Dyson 28 - thus showing the stochastic motion of the 


eigenvalues in our Bernoulli ensembles tends to the familiar motion of eigenvalues in 
the canonical Gaussian ensemble setting. 

The first example, in Section [5j concerns real-symmetric Bernoulli matrices. The 
second (Section [b]), concerns imaginary anti-symmetric Bernoulli matrices and the 
third (Section[7]) investigates the singular values of real non-Hermitian NxM Bernoulli 
matrices. These ensembles are closely related to the graph types ([!]) , (|m| and © 
respectively (see above). In the three cases we recover the JPDF corresponding to 
the fixed trace GOE, antisymmetric GUE (see |9,34 ), and real Wishart distributions 
(see e.g. [35] and references therein) respectively. Finally, in Section |8j we discuss the 
results and outline potential developments. 


2. Random walks on Bernoulli matrices 


Let Bn denote one of the unconstrained ensembles of N x N balanced matrices, as 
defined in the preceding section. Each matrix B £ Bn consists of entries By £ {±1}, 
of which djv of these variables are independent. Hence there are Bn = 2 djv matrices 
in the ensemble, with each one possessing equal probability to be selected. However, 
rather than simply selecting a matrix at random we employ an alternative method that 
allows us to recover this uniform distribution. The idea is to initiate a discrete random 
walk process over Bn- If this random walk is ergodic then a uniform distribution is 
recovered in the infinite time limit as the system approaches equilibrium. This walk 
is a discrete analogue of Dyson’s Brownian motion model for the Gaussian matrix 
ensembles. 

To start, let us suppose that each matrix represents a vertex in a large meta-graph 
(so called because it encompasses all graphs in our ensemble). Notice that each matrix 
is in one-to-one correspondence with the d/v-dimensional vector v(B) = (iq,..., i’d N ), 
Vi £ {±1}- Therefore the vertices of the meta-graph correspond to the corners of 
a djv-dimensional hypercube with 2 djv vertices. Two vertices are then said to be 
connected, or adjacent (denoted B ~ IT), if they differ by exactly one matrix element, 
or alternatively if their Hamming distance D{B,B') := i||u(B) — v(B ')||i is equal 
to 1. This connectivity relation thus means each matrix/vertex has d/v immediate 
neighbours. 

The analysis of random walks on the hypercube are well established and we refer 
the reader to [30 3ll, and references therein, for more details. Here, following [31], 


we define the transition probability that the walker either stays in the same place, or 
moves to any neighbour in one time step as equal to l/(djv + 1). Thus, if P t (B) is the 
probability to be at B at time t, then the probability distribution a single time step 
later is given by 

p t+1 {B) = [ e p t m = j±— Y, p ^ B 'y ( 2 ) 

N B'-.D(B,B')< 1 

This particular Markovian evolution ensures there is a unique stationary distribution 
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given by Poo(B) = 2 djv for all The authors of |3l| perform a comprehensive 
analysis of this approach to equilibrium by starting with a probability distribution 
P 0 (o) = 1 localised at the origin, which we choose to be v(o) = (—1,..., —1). They 
obtain an asymptotic expression for the total variation distance 


m = 


E 

B£B N :Pt(B)>2~ d N 


\P t (B)-2~ dN \ ~ Erf( e - 2c(t) /v / 8) 


( 3 ) 


in the limit of large d/v- Here Erf(a;) = f* dz e~ z /(27t) is the usual error function 
and eft) = (t — t C rit) /djy. It shows that V(t) remains very close to 1 until a critical 
time, calculated to be f cr it = jdN^ogdw, after which the probability density relaxes 
to the equilibrium distribution exponentially fast in the region t > f cr it and large d n 
limit. 

v(,, ~i; exp (“ 2 ^f)' (4> 

This exponential approach to equilibrium as t —» oo is controlled by the second largest 
eigenvalue (the largest taking the value 1) of the transition probability matrix g in 
©■ Due to the symmetry of the hypercube, the eigenvalues of g are equivalent, up to 
degeneracies, with those of the transition probability p, which describes the distance 
from the origin, outlined in Section [ 3 ] These are given by Aj = (1 — hJ+t) f° r 
j = 0, ... djv- Hence, we find the associated quantity C(t) = \\g t 5 0 —P oa \\ ~ (1— rfjV 2 +1 )* 
(where S 0 denotes the probability density that is equal to 1 at the origin and 0 
everywhere else) displays the same asymptotic approach to equilibrium In 

particular one observes this occurs at a rate controlled by the hypercube dimension 
djv and so we define a rescaled time 77 = t/dpf , which is the natural time-scale in this 
context. 


3. Induced random walks and entropic forces 


In Section [4] we analyse the joint probability distribution for the eigenvalues by 
studying how the random walk in Bn induces a corresponding random walk in the 
space of spectra. This projection amounts to a reduction of the information about 
the matrix. As a result, points in our ensemble, that were once distinct, become 
indistinguishable, or much closer, in the reduced space. This change of viewpoint 
creates effective forces that originate due to entropic, rather than mechanical, reasons. 
Consequently, one obtains a Fokker-Planck equation that describes how the random 
walker exhibits preferential concentration in certain regions in the space of spectra. 

Before we progress, as a prelude, we provide a brief example that illustrates this 
phenomenon: We compute the distribution of the Hamming distance of the random 
walker from the origin X(B) = D(B , o). This is in fact equal to Tr AA^ = ^I(A) 

(where A is the adjacency matrix). It is therefore the simplest instance in which the 
random walk relates directly to the spectrum. A numerical simulation of this walk is 
given in Fig. [lj where several random trajectories of X(B(t)), all starting at the origin 
o, reach equilibrium around some critical time f crit . Similar problems, in various other 
contexts, are well discussed in a number of articles [29-31,33], to which we refer the 
reader for more details. Here we shall outline the essential points. 


| In contrast, if one does not allow the random walker the chance to remain in the same place then 
there is no unique stationary distribution, since the hypercube is bipartite. 
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V 

Figure 1. Five simulated random walks over the space of matrices ( |21[ ), with 
A r = 50 and an initial matrix satisfying = 1 V i,j. Plotted is 

X(ij)/(d\- +1) as a function of 77 = 1/d ; y, with saturation to equilibrium occurring 
at approximately Ucrit = tcritMv = log (d N )/4 & 1.78. 


To begin, notice that X(B ) is given by the number of l’s in the code v(B), 
therefore the probability that the walker moves from B to another B' with X(B') = 
X1(B) — 1 is proportional to X. Similarly the likelihood of increasing X by 1 is 
proportional to dx — X. Given this transition probability is the same for all such 
matrices we drop the explicit reference to B and write 

X _d N -X 1 ... 

PX — >X — 1 j i i 7 PX —►Jf-f-l j i • Px — >X i i 

d N + 1 d N + 1 d N + 1 

with 0 otherwise. Therefore, if Pt(X) is the probability to find a vertex with a given 
Hamming distance X at time t, then the probability distribution is updated at the 
next time step via the following process 


Pt+i(X) - P t (X) = Pt(X - I)px- 1 -+X + Pt(X + l)px+i^x - P t { x)(l - px^x) 
d N 


2 (d/v + 1) 


(P t (X + 1) + P t (X - 1) - 2P t (X)) 


+ P t (X + l) 


X + l-d N /2 
dx T 1 


-P t (X-l) 


X-l- 


dx/2\ 


dx T 1 / 


( 6 ) 


The 


A detailed study of this discrete Markovian evolution can be found in 29 
spectrum of the time evolution operator (5|) is given by Aj = (1 — dN +]_ ), with 
j = 0,1,..., dx- As mentioned previously, tms spectrum coincides with the evolution 
operator for the uniform hypercube in Section [2j The authors of 31 therefore utilise 
the solution of (|6| in order to analyse the approach to equilibrium of ©• The 
stationary distribution for Pt(X) corresponds to the first eigenvalue of p (given by 
Aq = 1) and is reached when the cube is uniformly covered: 


Poo(X) 


1 

2 dN 


dx 

X 





N —>• oo. (7) 


The asymptotic estimate above is only valid for |A' — 4p| S o(d^) , and it gets worse 
when X is farther away from At the centre, the error is of order 0(d N 2 ). In the 
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continuous limit, as we approach equilibrium the probability density becomes normally 
distributed with mean p = /2 and variance cr 2 = d/v/4. It shows that the majority 

of our vertices on the hypercube lie in a region of size proportional to y/dN centred at 
a distance d/v/2 away from the origin. Thus our random walker is ‘attracted’ to this 
region simply because the number of vertices outweigh those near the origin. 

The result Q can also be obtained by taking the continuous limit of the equation 
(6) instead. This is done by introducing a scaled spatial variable £ = € 


[—\/djv/2, y/djy/2] and a scaled temporal variable p = t/dpf. Then, in the limit 
the discrete evolution equation approaches a Fokker-Planck (Smoluchowski) equation 
describing an Ornstein-Uhlenbeck process 



dP&V) _ l &P&T,) d(ZP(Z,T,)) 

dp 2 df 2 <9£ 


Setting the LHS equal to 0 and solving for P(£, oo), one obtains 0 in its scaled form. 
Several comments are in order at this point: 


(i) Choosing a different scaling of X, such as x = X/dx — 1/2 with x G [—1/2,1/2] 
might have appeared natural in some instances, and would lead to P(x, oo) = 
6(x), emphasising the concentration of the probability to the centre of the X 
interval. However, one then loses the important information that in the limit, the 
variable £ behaves as a Brownian process. 

(ii) The scaling of X with y/dN is the natural scaling in the present context. It is 
equivalent to scaling Tr AA^ by N and hence A by y/N. For this reason we employ 
the same scaling when investigating the spectra of Bn in Section [4] Similarly, the 
time scaling p = t/dn will also be used in Section |4j 

(iii) The second term in the RHS of the Fokker-Planck ([8]) is usually interpreted as due 
to an effective drift force, which takes here the form F(£) = —2£ as if it derives 
from a confining harmonic potential. Clearly, in the present context it arises 
purely because of the combinatorial (entropic) preferences in the hypercube, which 
create the illusion of a confining harmonic potential. This mean drift towards the 
centre £ = 0 will also appear when we investigate the stochastic nature of the 
eigenvalues with respect to perturbations of the matrix B. Note also that equation 
Q plays a crucial role in Dyson’s Coulomb gas model for random matrices 
which he uses in order to introduce the assumed Brownian motion of matrix 
elements for the Gaussian ensembles. 


28 


4. Spectral statistics 

We now turn our attention to the induced random walk in the space of spectra. The 
formalism will be presented in general and then applied systematically to specific 
ensembles of Bernoulli matrices in Sections [5j [6] and [7] Following from the discussion 
of point in the previous section we introduce the scaled matrix B = B/y/N, to 
which we associate an ./V-dimensional vector of eigenvalues (or singular values as will 
be the case in Section [7]) 

cr(B) = (A !(B) < A 2 (B) <...< A N (B)). (9) 

This ensures the eigenvalues A„ are 0(1) (and thus the mean level spacing is 0(1V _1 )). 
Furthermore, in each ensemble, we shall find that every matrix B G Bn satisfies 
Tr(H^ B) = R 2 n , where Rn is an iV-dependent constant (to be calculated) of size 
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0(N 1 / 2 ). Hence, the set of spectra, which we denote by Ejv, forms a discrete and 
finite subset of the continuous space Sn = {s € : s„ < s„+i, ||s|| 2 = Rn}- 

Suppose now we initiate our random walk, starting at and moving through 

the meta-graph according to the transition probability ([2]) . What are the dynamics of 
<t(BW) g Ejv as we follow this walk over the integer-valued time t? Fig. kl| provides 
a numerical simulation of a ‘typical’ realisation of this walk: Every eigenvalue in u 
is plotted as a function of the scaled time r\ = t/djy. The main point to observe 
here is the clear separation between the microscopic, mesoscopic and macroscopic 
times-scales. In the main figure we see how the eigenvalues approach equilibrium on 
scales of order of the system size t = O(dN) (in contrast to Fig. [l] in which it was 
dw log(djv)/4), however before reaching this stage coherent trajectories already emerge 
due to the drift forces associated with the ‘repulsion’ of eigenvalues. This happens on 
a time scale t = O(N), which is more evident in the inset picture. In the microscopic 
regime we have fluctuations arising from a change in the spectrum Scr due to a single 
time step, which we denote 5t = 1. Between these sits the mesoscopic regime which 
characterises the point at which these spectral trajectories start to emerge. Here 
we obtain fluctuations Act corresponding to small (but still larger than microscopic) 
changes in time At = 0{N C ), 0 < c < 1. In the remainder of the paper we shall 
continue with this convention of using S to denote variations due to a single step and 
A to denote variations due to At time steps. The clear hierarchy of time-scales is one 
of the key ingredients which enables the construction of a Fokker-Planck approach to 
describe the spectral evolution which emerges due to the underlying random walk. 



o 0.5 l , 1.5 2 

T >=dk 


Figure 2. A single instance of the spectral evolution during the same random 
walk process as Fig. ^ Plotted is the entire spectrum as a function of rj. The 
inset provides a close up view of the same process. 

The evolution of the spectral distribution on the microscopic time scale starts 
with a Markovian master equation analogous to Let 7 r(cr,£), be the probability 
to be at some point a G Ejv at time £, subject to some prescribed initial conditions[§] 

§ Note that at this level of the description, both the time t and the space coordinates <r are discrete 
variables. 
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Then for b(cr) = {B £ Bn '■ &{B) = er} we have 

, 1 1 

7r(cr,f + 1) = T r7 - T T ^ 


|6(cr)| 

1 V n Bdb(c 


dN + 1 


n{cr(B'),t), (10) 

B':D(B,B')<1 

oo, 7r(cr,f) will converge to the 


where D(B,B r ) = D(B,B'). In the limit t 
stationary distribution associated to the JPDF for the eigenvalues of Bn- However, 
in contrast to the example in Section [3] there is no closed form expression for n(cr,t), 
or even 7 r(<x, oo). Instead, the best we can hope for is a convergence as N —> oo to a 
limiting distribution, analogous to the Ornstein-Uhlenbeck process ^ in the previous 
section. Unfortunately, the obvious relationship between ([8]) and its discrete process 
([6]) has no analogue in the context of the induced random walk in the space of spectra. 
For this reason, we introduce an alternative derivation that overcomes these failures. 
It is an adaptation of Kolmogorov’s method 32], which is eloquently presented in 
Section 8 (c) of Wang and Uhlenbeck (33 (other approaches using essentially the 
same arguments are also outlined in 36 37 ). 


In the classic example of a random walker on a ID line (see e.g. 37 Section 3.8.2), 


the Brownian motion formalism emerges from its underlying microscopic dynamics 
because the motion of the walker is described exactly by a sum of independent random 
variables. In other words the central limit theorem (CLT) guarantees for a time 
window of At, the transition probability is Gaussian, up to an error (At) -1 / 2 due to 
the Berry-Esseen theorem (see for example 17:; 38 for various versions of this theorem 
and previous applications to RMT). Whilst we shall not invoke the CLT explicitly 
in our context, it is helpful to keep this line of thought in mind, as it underlies the 


transition from the microscopic discrete view provided by (10) to a coarse grained 


Fokker-Planck description. This is made possible due to the clear separation between 
the microscopic and macroscopic time scales. On the intermediate (mesoscopic) time 
scale At ~ N c , 0 < c < 1, our evolution will be given as an accumulation of, 
approximately, independent random steps. At the same time, the spatial resolution 
is deteriorated due to the increasing width of the distribution acquired during (At), 
which will be shown to be of order N ’A . Thus, the temporal and spatial coarse 
graining have opposite dependence on At, and an optimum description should be 
sought by the freedom of choice of the parameter c. The embedding of the discrete 
evolution in continuous space-time is made possible by the fact that the macroscopic 
time scale is O(djv) = 0(N 2 ), so that on this scale At can be approximated by 
a differential. However, one should bear in mind that as long as N is finite, the 
continuous space-time evolution is not precise and the errors accumulated by various 
approximations will be discussed in detail. 


4-1. The Fokker-Planck equation 

Evolving our random walk by a single time-step equates to choosing an entry B pq at 
random from the matrix B and changing its sign. This produces a new matrix B' , with 
the difference matrix SB = B' — B of either rank 1 or 2. The form of SB with respect 
to this change is given, for example, in (221. In At time steps the random walker 

a maximum Hamming distance 
,ua t) i n our hypercube 


is capable of moving from B to another matrix B ', 

D{B, B') = At away. This corresponds to a path uj = (wi, 

Bn, where each denotes the swapping of a particular matrix element pgjjj] The 


|| There is also the possibility of remaining in the same place, which is treated more thoroughly in 
Section E3 










A random walk approach for Bernoulli ensembles 


10 


perturbation associated with u> is simply the addition of all single step perturbations 
along this path 

At 

AB U = B' — B = ^ 5B Ur . (11) 

r—1 

The probability to move from the matrix B to B' = B + AB in At time steps 
is denoted pAt(B —>• B') and given by g At , as defined in ([ 2 ]). Employing this coarse- 
graining in time the microscopic Master equation ( |10| becomes (abusing notation 
slightly and writing in terms of the rescaled time variable rf) 

7r(cr, 77 + A77) = ^2 PA v (c' a)n(a',r]). (12) 

ct'GEn 


Here the transition probability in the space of spectra Ejv is given in terms of the 
underlying evolution in the space of matrices, i.e. 

PAr,{cr -*• er') := —^ ^ p At (BB') e r'), (13) 

1 ' '' flefe(er) B':D(B,B')<At 


where 5(a,b) is the Kronecker delta. However, in Section 5.1 we show (in the large 
N limit) this transition is independent of the initial matrix B and so rather than 
summing over all matrices in b(cr) we can choose any particular matrix B that satisfies 
er = cr(B). 


PAriicr -t cr') := ^ p A t(B -t B') l5(e^(H , ), er'), (14) 

B':D(B,B')<At 

The next step in the derivation requires the construction of the time-averaged 
moments of the evolution, 




E(Ilti 

Arj 


— A 77 ^ ^i)P&v( cr 


E 

'es 


(15) 


i =1 


had 


from ([14 1. 

If we 
ensembles 

moments survive the limit A 77 —> 0 33 


28 


a continuous process, as in Dyson’s prescription of the Gaussian 
the derivation would proceed by showing that only the first two 

In the present case our process is discrete 


and therefore for any finite N there is bound on how small we can take A?y. For this 
reason, one can only show that for some finite N the first two moments are (to leading 
order) independent of Arj and that higher moments will depend on higher powers of 
Ap. This persistent presence of the higher terms is the main source of error in the 
transition from the discrete to the continuous description. 


Following the approach in 33 , we introduce a test function h(cr), which vanishes 
on the boundary Sn, together with its derivative. The requirements on the smoothness 
of h(cr) and on its domain will be specified later. We now study the coarse grained 
time evolution of the expectation value 


K~ E ft ( <T ) 7r ( <T > r ?)- 


<tGEj\ 


Following (12 1 we have 

h v +A v = ^2 /l ( <T ) 7r ( cr ; P + A 77 ) 

trGSiv 


(16) 


^2 Mr/( cr -t (t')'k(ct , rj)h{a'). 
<T,<r'GSiv 
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Expanding h{a) as a Taylor series, dividing through by Arj and performing the sum 
over cr' we obtain, using (15), 


Air 


Y h ^Ar, = Y 7r ( <T >»?) 






v-^ ,, dh 


+Y Mvi * 


d 2 h 


2 dX v d\„ 


(17) 


where Air = i r(cr, rj + Arj) — 7r(cr, rf). 

It was shown above that the coarse grained temporal description induces naturally 
a corresponding spatial coarse grained representation of ir{a,r i /), and therefore it is 
advantageous to convert the sum over the discrete set of spectral points E^r to an 
integral over the covering domain Sn. In contrast to the previous section, in which 
the points X were arranged in a lattice, the spectra cr are distributed through Sn in 
some random fashion. We therefore assume the existence of a density function q{a), 
so the set Ejv effectively appears to be a sample of |Ejv| points chosen randomly from 
the distribution q{a). Given this assumption the expectation of some function / with 
respect to q{a ) is 

[ da f{a)q{a) « — l - r f{a) = E Sjv (/). 

o-e£j\ 


' S N 


-‘N] 


This approximation is the premise of the Monte-Carlo integration technique. If indeed 
the points a can be assumed to be random distributed points according to q{a) then 
the error is essentially given by the law of large numbers, i.e. yWar (/)/|Ejv|, where 

Var {f) = jY- J2 (/(<t)-E S jv (/)) 2 

crGS 


->N 


N 


Importantly, this error does not depend on the dimension of the space, only on 
the number of points |Ejv|- Asymptotically this is |Ejv| ~ 2 dN /N\, which grows 
exponentially as N increases, meaning the error decreases very fast. Moreover, the 
density of spectra is 

|Ejv| 2 n ~ 1 T{N/2) 2 jv ( jv + 1 )/ 2 


„JV 2 In 2—zN In N 


Z > 0. 


\S N \ 2ir N / 2 (N + l) T(N + l) 

This leads to an exponential increase in the number of spectra per unit (hyper) volume, 
allowing us to take functions /(cr) on a scale far smaller than the typical spacing 
between eigenvalues. We therefore write 


' S N 


da f(a)ir(a, rj)q(a) 


cr e£j> 


7r (°' 1 ?7 )f(a), 


with ir(a, rf) a smooth approximation of || tt(ct, rf). In fact the main source of error 
will come from the higher terms in the expansion and this sets a much greater 
restriction on our test functions. Writing Q(a,rf) = it (a ,rf)q( y a) and replacing the 
sum in (|17| by an integral we arrive at the following expression 


/ dah{a)^ 

ls N A v 


' Sn 


da Q(a,rj) 


v—> Oh 1 v—^ 

E M -gr + 5 E M ’ 


Iv'/i, 


d 2 h 
v>x d\ v d\ u 


+£{N,h{a)).{18) 


The error £(N,h(a)) depends on a number of factors and will be discussed below in 
detail. If we replace the time difference in the LHS side of (181 by a derivative (by 
sacrificing an error of 0{Arf) = 0{N C ~ 2 )), and integrate by parts the RHS of (18) 
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(where we make use of the requirement that h and its gradient vanish on the boundary 
of Sjv) we get 



dcr h(cr) 



d{M„Q) 

<9A„ 


1 d 2 (M Vfl Q) 
2^ dXudXy. 

Vfl ^ 


£(N, h(cr)). 


(19) 


Consider the above relation where the error term is replaced by 0. In order to satisfy 
this equation independently of the actual choice of the observable h, the expression 
in the curly brackets must vanish. Thus Q(cr,p) should satisfied the Fokker-Planck 
equation: 


dQ _ d(M v Q) 1 ^ <9 2 (M l//I Q) 
dp A-j g\ v 2 d\ v d\n 


However £(N, h(cr)) does not vanish. At best, under certain smoothness constrains 
on h the error term vanishes as a negative power of N in the limit N —> oo. Hence, for 
finite N the evolution of the spectra under the random walk is described by the solution 
of the Fokker-Planck equation within an error arising from a number of factors, as 
detailed below. We also note that it will transpire that the off-diagonal second order 
moments M vfJj (with v ^ p) provide a lower order contribution than their diagonal 
counterparts M vv . Therefore, it is amenable to remove these terms, and place their 
contribution within the error term £(N, h(cr)). Doing so, we shall obtain a Fokker- 
Planck equation in Sections[5j[6]and[7]with only a single sum over v, which is consistent 
with the resulting equation in Dyson’s approach [281. 


).2. Error estimates 


There are several sources of error which occur in the transition from the discrete to 
the continuous description of the spectral evolution for finite N: 


(i) The replacement of discrete sums (over space) and differences (in time) by 
integrals and differentials, respectively. The error induced by the former was 
discussed above and shown to be exponentially small, while the latter is of order 

0( ^ n -{2-c)) _ 


(ii) Corrections to the leading order expressions for the moments M v and M vv . In all 
three ensembles we consider, these corrections will arise from higher terms in the 
perturbation expansion and will require assumptions regarding the delocalisation 
of eigenvectors and the rigidity of eigenvalues. The conditions under which the 
assumptions are valid have recently been stated (see e.g. fl7[[l9p9| and references 
therein for details) and are known to be satisfied with a high probability (tending 
to 1 in the limit of large N). 


(iii) Off-diagonal second order moments M„„ (with v p), as described above. 

(iv) Moments in of order 3 and higher. 


This last source of error will in fact provide the greatest contribution to £ ( N , h{ar)) 
which, as seen in © is coupled to higher derivatives of h(cr). To ensure convergence 
(i.e. for £(N, h(cr)) to decrease in N) we must therefore restrict our observables h(t r) 
to a certain class with appropriate resolution determined by the following restrictions: 

1. The support of its Fourier transform is restricted to a ball of radius k = 
0(N l+a ), a > 0, which means the derivatives of h , in any direction, are of order 
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fo( fc ) = 0(N k ( 1+a ^)h. In particular, this choice ensures that the resolution is always 
better than the typical spacing between eigenvalues, which is 0(N~ 1 ) 1 as the width 
of h is of the order of k~ 4 . 

2. The test function h.(cr) is defined over a spectral set T = {A yi }" =1 which 
consists of n eigenvalues with n = 7 > 0. This limits the size of the spectral 

interval under study to be smaller than the entire spectral support but to include a 
large number of eigenvalues. 

Under the restrictions above, it is possible to estimate the magnitude of the 
neglected terms. For the purposes of clarity we defer the computations to Sections [ 5 TT| 
and Appendix A Here, we simply mention they are consistent with the results one 


would suspect from using the CLT in combination with perturbation theory. Namely, 
that all moments beyond first order are consistent with a Gaussian kernel pAri(cr, &') 
comprising independent AA„, each with mean E(AA„) ss E((z/|<5-B|i/))Af = 0(1)A?/ 
and variance E(AA 2 ) ss Fi{{v\5B\u) 2 )At = 0(N~ 1 )Ap (the exact expressions are 
given in Section |5.1| ) . Only the first moment will deviate from this prescription due 
to the influence of the second term in the perturbation formula (231, though it will 
remain of the same order. Therefore, summarising the contributions we have. 


First order : The first moment is of order 0(1), so 

[dcrQY^M,^ = [ dcrhQ[0{N 1+ ' 1+c 




using that h' = kh = 0(N 1+a )h 

• Second order: The main contribution to second order comes from the diagonal 
moments M = 0(N~ 1 ). Therefore 


d 2 h 


= I dahQ[0(N 2+2, 


: q +7—1 


uer 


)], 


using that h" = k 2 h = 0(N 2+2a )h. Note the off-diagonal terms will give an 
overall contribution of 

J dcrQ Y = / d <r h QP( N2+2a+2l+C -*)l 


which is much less than the diagonal terms. 

Third order: The main third order term comes from moments of the form M, 


Therefore we take 


VI/fl’ 


d 3 h 


daQ Yj 


= J dcrhQ[0(N c ~ 3+3+3a+2 ' 1 )] 


• Fourth order: The fourth order term comes from moments of the form M, 


vv mi 


d A h 

QXWXl 

v u 


dcrQ Yj 


= J dcrhQ[0(N c ~ 4+4+4a+2 ' 1 )\ 


We require the third order to be less than the second. Hence l + 2a + 7 > c + 3a + 27 , 
or a + 7 + c < 1. Similarly, for fourth order we get l + 2a + 7 > c + 4a + 2y, 
or 2a + 7 + c < 1. The latter provides a slightly tighter restriction on a. Thus, if 
2a + 7 + c= l — e, the error estimated for ( |l6|) d ecreases as N~ £ . 

Ideally, we would like to translate the result (|19| into a statement comparing the 
JPDF for our Bernoulli ensembles, with the resulting stationary distribution Q(cr , 00 ) 
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of the Fokker-Planck equation (20). For instance, a weak convergence estimate given 


by the difference in the expectation of our observable. At present we do not know 
of any way how to achieve this, however we believe this difference to be of the order 
of | M <T ) 7r ( <T ; °°) — Ss N c ^ T M cr ) ( 3( c U °o)I = 0(N~ e ), which comes from the 

ratio of the second and fourth moment contributions. This is because after setting 
the time derivative to 0 and using the estimates of the derivatives of h, we obtain 
a perturbative form of the Fokker-Planck equation, i.e. LppQ = xQ , where Lpp 


denotes the operator on the RHS of (20) and x = 0(N~ e ). A precise error estimate 
will necessitate a careful stability analysis of the Fokker-Planck equation, which is 
beyond the scope of the present article. 


5. Real symmetric Bernoulli ensemble 

Our first example consists of the following matrix ensemble Bn, in which the matrices 
B e Bn have elements 

±y/l/N p < q 


B pq — 


B q p p> q 


( 21 ) 

± s/2 /N q = p, 

chosen with equal probability. The addition of the diagonal elements in comparison 
to the balanced matrices 0 serves to simplify the proceeding calculations but is not 
crucial for the result. It does mean the number of free parameters, and therefore the 
dimension of the as soci ated hypercube is given by dN = N(N + l)/2. Moreover, as 
one can verify that Tr(R 2 ) = 2 dw/N = 0(N) for all B £ Bn- 


remarked in Section 


4.1 


The ensemble (21) falls into the class of Wigner matrices and hence one can show that 
in the limit of large N, the mean spectral density of the ensemble approaches that of 
the Wigner semicircle p( A) = \J 4 — A 2 /(27r) (see e.g. 41 for a definition of Wigner 
matrices and proof of this fact). 


5.1. Evaluation of moments 

Let us suppose that we have two matrices B' ~ B that are neighbours in our meta¬ 
graph. B' is thus obtained from B by switching a single element or, equivalently, by 
adding a rank 2 (or rank 1 if it is a diagonal element) matrix SB to B , given by 

~2B pq (\p)(q\ + \q}(p\) p < q 
-2B pq \p){p\ P = q- 

Here we use the ket notation | p) to denote the column vector with a 1 in the p- 
th component and 0 everywhere else, and the bra (p | denotes the transpose of |p). 
Importantly this perturbation depends on the original matrix B, and therefore so will 
A B in (11), provided that our random walker does not move back and forth in the 


SB (pq) = B' — B = 


( 22 ) 


same direction of the hypercube. Or, in other words, all uy belonging to the path 
to = (wi,..., uiAt) are different. However, this situation occurs with probability 1 in 


the large N limit (see Subsection 5.1.1), since the number of possible choices dN, far 


outweighs the number of steps of the random walker At. 


The perturbation (11) in turn induces a change in each of the eigenvalues AA y . 


Unfortunately this change cannot be calculated exactly but it can be approximated 
using the well-known perturbation formula 

AAj, = (i/\AB\u) + S/ + higher order terms , (23) 

z —' A,, — A,, 
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with the third order term given in (A.l). In complete analogy with Dyson’s Brownian 
motion, our aim is to use this expansion in order to find expressions for the moments 
(5.1), which will come from the linear and quadratic terms in the matrix perturbation 
A B. The difference is that we are dealing with a discrete process and therefore we 
must be careful regarding higher terms in the expansion, whereas in Dyson’s case the 
higher terms are justifiably neglected as the time differential can be made arbitrarily 
small. Nevertheless, (111 gives us an indication why this might be viable. Looking 
closely we see that, even though the elements of the matrix A B will only take one of 
two values, each of the individual random perturbations SB are (almost) independent 
and identically distributed. Therefore, by reducing to the space of spectra, terms such 
as (i/\AB\i>) will become a sum of iid random variables and hence, due to the CLT will 
approach a Gaussian distribution with appropriate mean and a variance proportional 
to At, which is the essence of Brownian motion (see e.g. 1(33,, 36]). 

This observation also imposes a restriction on the range of At for which the 
perturbation formula provides a good approximation to A\„. Since (v\AB\v) is the 
main term in this approximation, we require this to be much less than the typical 
spacing between eigenvalues, given by A = 4 /N . Therefore, since the variance grows 
linearly in At the standard deviation is 

\/Var(i/| AB\v) 

~A 



where we have used that E((i/| AB\v) 2 ) = 8/(Ndx) (see Subsection 5.1.2) and 
d/v ~ N 2 /2. Thus we require At <C N, which is comfortably satisfied if we take 
At = 0(N C ) with c < 1. We shall show below this same restriction also emerges from 
different reasons. 

In the following we shall calculate the expectation of terms both linear and 
quadratic in the matrix perturbation A B, for example (v\AB\v) and {v\AB\p) 2 . 
These will be calculated explicitly and will provide support for the line of thought 
discussed in the introducing remarks to section [4] - concerning the relation to the 
CLT. Collating these will allow us to obtain expressions for the first two moments M v 
and M vfl (see (15)) which in turn allow us to obtain a Fokker-Plank equation and 
therefore the stationary JPDF. This is presented in Subsection |5.2[ The estim ation of 
higher orders, specifically cubic and quartic terms, are left until Appendix A[ 


5.1.1. Linear terms We begin by calculating the expectation of the first term 
(v\AB\v) in (23) due to a random walk of At time steps. Given the linearity of 
this expression we can directly calculate the expectation in the matrix perturbation, 
E(A B), which, using (141, is given by 

^{AB) = Y j Pm{B^B')AB{B')=Y j I ^1 Y. ( 24 ) 

B' x=0 A ' B':D(B,B')=X 


Here P&t(X) is the transition probability to go from any matrix B, to another B' a 
distance X away, as described in Section [3] The probability to reach all such matrices 
is the same and thus we simply divide by their number Afx- Alternatively, this can 
be viewed as summing over the set fix of all ordered paths of length X. These are 
paths u> = (wi,... ,u>x) in which all uj r are different and arranged in, say, decreasing 
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order and each such co £ fix corresponds to a unique B ', a distance X away from B. 
Hence \Nx = |Hx| and 


E(AH) = 


^ PAt(X) 


At 


I fix I 

x=o 1 1 


E AS " = £ 


uj^Clx 




PAt(X) 

|fi 




■X 


oj^Qx r=l 


E^ v ’( 25 ) 


where we have inserted the summation over single step perturbations © and, in 
comparison with the stationary distribution 0, 


ia Y i = 


djv 


(26) 


The double sum in (25) can be written as a single sum over the nearest neighbours 

of 5, i.e. E r = 1 SB air = J2 p < a ^x U±J ' '> wirac 


, p<q 5B pq , where d>^E denotes the number 


of paths ui £ fix that contain the element pq. However, this must be the same for 
all (pq), so without loss of generality we can write = d>E which is given by 

conditioning (pq) to belong to w and choosing the remaining X — 1 elements from the 
dx — 1 possible options, i.e. 




d/v — 1 
X - 1 


X 

dx 


dx 

X 


= x— i a* i 

dx 


Therefore removing d>^ outside the sum we obtain 


At 


E(A.B) = E 


x=o 


PAt(X) 

l n x| 


V SB pq = 


"x zA 

P<Q 


1 

dx 


At 


E XP M(X) 


x=o 


/A 

P<9 


S B pq 


(27) 


(28) 


We are left with two seperate contributions to determine; one from the expected 
distance of the random walker away from the initial matrix B after At time steps, the 
other from the aforementioned nearest neighbour sum. Let us start with the former. 
We know this cannot exceed At, therefore we can write 


At 

AtP At (At) < E XPAt(X) < At 

X=0 

for all At > 0. It thus suffices to estimate P A t(At) for At <C dx- From Section[2]we 
know that at each time step the random walker has (dx + 1) options to either stay in 
the same place, or move to an immediate neighbour. Thus after At steps there is a 
total of (dx + 1) A * possible paths the random walker can take, with Af!|f2At] of these 
options leading the random walker to a matrix B' that is a distance At away from B. 


Hence, for At <C dx, we have P At (At) - 


(dx + l) At dx — ^ ~*~ 

j=o 

= 1 + 0(At 2 /d N ) = 1 + 0(N 2c ~ 2 ) and so 


(29) 


E(AH) 



This leaves us with the second contribution from the nearest neighbour sum. 
Using the expression for the single step perturbation (f22j) we have 


E SB™ = -2 E B pq (\p) (q\ + \q)( p \) - 2 E B pp \p)(p\ = —2B.(30) 

p<q p<q p 
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Therefore 


E(HA^)) 

Arj 


= —2{v\B\v){\ + 0{N Zc ~ z )) = -2\ v + 0(N Zc ~ z ), (31) 


where Ary = A t/d^ and A„ = 0(1). Up to a correction, this result is exactly what 
one would expect for the Gaussian ensembles, showing there is an overall ‘force’ on 
the spectrum that acts to bring the eigenvalues closer to 0. However in contrast 


to Dyson’s case 28 it arises not because of an imposed Gaussian distribution of the 


matrix variations, but because of the entropic forces experienced by the random walker 
as it moves within the hypercube, just as in Sectio n [3] The restriction of c to the 
interval 0 < c < 1 is necessary for the error term in (31) to tend to 0 when N —> oo. 


5.1.2. Quadratic terms We now turn our attention to those instances in which 
we have contributions quadratic in AH. There are two such scenarios. The first 


occurs in the second term of the perturbation formula (23), the second in the leading 
contribution to M^„. The former concerns quantities E(|(^|AH|p)| 2 ) with v /x. 


whereas the latter those of the type ~E((u\AB\i/)(fi\AB\fj,)). We start by evaluating 
the first case in the same manner as for the linear terms. The expectation is therefore 


given, analogously to (25), as 


At 


E(|HAH|/i)| 2 )= £ 


X=0 


PAt(X) 
|fl 


A' 


A' 

E E 

cuGQx r,s =1 


y\6B“r\pi){n\5B“°\v) 


where, we have used the expression |llj) for the multi-step perturbation of A B. Here 
we must now split the double sum in terms of diagonal and off-diagonal parts, i.e. 
J2 rs = Er+Er^s- The reason being that the diagonal term, as before, consists of 
counting those ordered paths u> that contain a particular element ( pq ), whereas in the 
off-diagonal term we must count the number that contain two distinctive elements. In 


fact the expression (271 generalises quite straightforwardly to determine the number 


of ordered paths of length X with l distinct elements, to 
■At) . I d.\ I 




x-i 




2 — 0 


(32) 


Using this we find 
E(|HAH| M )| 2 ) = 


^ PAt(X) 
x=o 


I^A’I 


$ 


Eiw<^»r 


d> 


xE' 

e^e' 


u\SB e \p)(p\6B e \v) 


,(33) 


where e denotes one of the djq independent elements pq of the matrix. The diagonal 
term can be evaluated in the same way as before, to get 


1 


d AT 

JV A=0 


At 

x " XP At (X)^ 


p<q 


{v\SB pq \ia)\ 2 = 


At v > 

dN ■“ 
p<q 


\{v\8B pq \n)\ z (\ + 0(N Zc ~ z )). 


(34) 


Summing over the nearest neighbours then gives 


p<q 




= EiA> 

p<q 


(iy\p){q\p) + B qp (iy\q)(p\p)\ 2 + J2 B 1 P \W\P)(P\^)\ 


4 

N 


E 

.p<q 


{ V pPq + 2 VpPqVqPp + VqP-p) + ^ ' 2 


— pj EK^ v pt L q v qt L p) ~ 


(35) 
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where we have used that v p. The off-diagonal term gives the following correction 


At (At — 1) 
d N (dN — 1) 


Eh 5B(pq) \d) 

2 

-Eih<^V)| 2 

p<g 

p<q 


= O(Ar) 2 


N’ 


where again, we have used that v ^ p and the prefactor comes from 
Exl 0 PAt(X)^/\i}x\ ~ -PAt(At)d>^/|f2At|- However this correction is less than 
the one arising in the diagonal term (p4|), thus 


E(IHAH|/i)| 2 ) 

Arj 


±(l + 0(N 2c - 2 )) = ^+0(N 2c ~ 3 ). 


(36) 


The second scenario involves terms of the form E((i/|AS|^)(/i|A£>|/x)). Therefore, 


starting once more from (25) 


At 


E(( 


v\AB\u)(p\AB\p)) = E E E 

x=o ' x ' Luen x r,s= i 


x 


\v\5B^\v)(pi\5B^\p) 


and proceeding as in the former case we get a diagonal term 

Ar,Y(^B pq W)(ti\SB^\p}(l + 0(N 2c - 2 )), 

p<g 

where it can be shown that 

Similarly, going through the steps above the off-diagonal part is given by 


0(Ar, 2 ) 

_P<9 

which equals 


[ E(^I^ (P9 V) E<^l ZB^\p) 

p<g ) \p<q 




0( Arj 2 ) 




= Ar 1 0(N c ~ 3 ) 1 

s + 0 (N 2c ~ 3 ) v = h 


(37) 


using A„ = 0(1). Therefore 

E((u\AB\v)(p,\AB\n}) _ w 

AT] \ 0(N C ~ 3 ) 

The results for the linear and quadratic terms from Subsections |5.1.l| and |5.1~2| 
together with calculations in |Appendix A.l] and |Appendix A.2[ provide expressions for 
the first two moments, together with error estimates (which rely on results concerning 
the delocalisation of eigenvectors (see e.g. 17 20 21 and references therein). In 
particular, for the first order moments we find 


Mv ~ ~ 2A " + ]y E 




A„ — A, 


+ 0(N c+q ~ 1 ). 


(38) 


Here the correction comes from the third order term in the perturbation expansion 
(A.l). The numerator of which is estimated in Appendix A.l whilst the denominator 


we estimate using the spectral rigidity properties of Wigner matrices provided in 17 
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This says that the spacing between eigenvalues is 0{N~ 1 ~ q ) with high probability, 
where q is some small positive constant. Similarly, the main correction to the moments 
M vv arises from expanding out the perturbation formula (231 and evaluating the cubic 
terms in the following 

y- E((v\AB\i')(i'\AB\fj,)(fi\AB\v)) 

^ A, - A„ 


Again we may use results in | Appendix A.l[ which show the numerator is of the order 
of 0(N C ~ 3 ), bringing us to 


M vv = - + 0(N c+q ~ 2 ). 


(39) 


5.2. Stationary solution 


If we take the leading order expressions for the moments (381 and (391 and insert 


these into (201, we obtain the same Fokker-Planck equation as Dyson 128 , which 


describes the stochastic motion of the eigenvalues in the GOE ensemble. However in 
our case we must condition the eigenvalues to obey a fixed trace, since Tr(H 2 ) = d N /N. 
Hence, setting the time-derivative to zero and solving for the stationary distribution 
Q(cr) = Q(<x,oo) we recover the JPDF for the fixed trace Gaussian ensemble |j9 


Q(rr) = C N H |A„-A M | <5 


/XO 


; Xl - 2d N /Nj , 


(40) 


where Gat is an appropriate normalisation constant. The properties of this ensemble 
are studied in various articles, where it is proven that the mean spectral density 
approaches the semi-circle law, and that the deviations for finite N follow the Tracey- 


Widom theory 34 . 


6. Imaginary anti-symmetric Bernoulli ensemble 

In our second example we investigate Hermitian Bernoulli matrices H with purely 
imaginary elements. If one takes the adjacency matrix of the tournaments graphs 
outlined in point jiu| in the introduction then we have a direct relation via the 
transformation H = iB = i(2A — J + I)/sqrtN. 

±i/\/N p < q 


H pq — 


II. 

0 


qp 


p> q 

q = p ■ 


(41) 


We would now like to discern the stochastic motion in the eigenvalues of the ensemble 
of H and from this extract the JPDF for the eigenvalues. In this instance all 
the matrices satisfy the anti-commutation relation H* = —H, which implies the 
eigenvalues come in pairs {±A„} with associated eigenvectors \v) and \v*) that are 
complex conjugates of each other. This means if N is odd there is a single eigenvalue 
denoted by Ao equal to 0. Note also the slight difference in this scenario that, due to 
the lack of diagonal elements, d/v = N(N — l)/2. 


We would also like to mention that recently, it has been proven in 48 , using a 


different approach based upon the work in [ 13f}16] , that all order correlation functions 
for the eigenvalues of H converge weakly to that of the associated Gaussian ensemble 


JPDF of (44). 
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6.1. Evaluation of Moments 

The evaluation of the moments follows in exactly the same manner as in Section |5.1| 
and therefore many details will be omitted. Using the estimate of the random walker 


a distance At away from the origin (|29| we once again find for the linear terms in AH 
E(A H) = (1 


(D(Ar]))A v y"SH p<1 . 


p<q 


Here the change in the matrix H induced by switching the pq -th element is SH pq = 
—2H pq [\p)(q\ — |q)(p|] is slightly different from the example in Section 


21 


However 


after summing over all elements we recover an analogous expression to (30) 

' SH P q = -2j2 H Pq [\p)(q\ - l?)(p|] = -2 H. 

p<q 


E< 

p<q 


Therefore, 


E(HA H\v)) 
Ap 


= -2A „ + 0{N 2c ~ 2 ). 


For the quadratic terms we recover an analogous expansion (331, with both 


of the form 


diagonal (341 and off-diagonal terms. The former of which gives a leading contribution 


E(|(^|A1T|/x}| 2 ) = A? 7 (l + 0(N 2c ~ 2 )) ^ \(v\SH pq \p}\ 2 , 


p<q 


with 

Ei 

p<q 


[v\5H Pq \p)\ 2 = — Bl^| 2 K| 2 - VpVqHpHq) = 


4 /N v ± p* 
0 v = p* 


Here we have used that \H pq \ 2 

,*2 


= 1/N and since \v) and \v*) are orthogonal 


eigenvectors J2 P = 'fZp u p = 0- Therefore, using the perturbation formula (231, we 
have, very similar to the previous section, 


Mv ~ M E 


N 


+ 0{N c+q ~ 1 ) 


and 


M v 


= — +0(N c+q ~ 2 ). 


(42) 


(43) 


The errors once again appear from higher terms in the perturbation formula, which 
require knowledge about eigenvector delocalisation and eigenvalue rigidity. 


If again we take the leading order expressions in (42) and (431 and insert these 


into the Fokker-Planck equation (19) we can obtain an expression for the JPDF 


Q{cr) = Q(er,oo). In contrast to the previous section, the eigenvalues come in pairs 
{±A„}, which we group together. This means if N is odd then we have the labels 
v = 1 ,..., (N — l)/2 so the expression for the first moment can be written as 


M v - —2A„ + E 

L 


1 


1 


Aw A,. 


Aw + Xp 


N A„’ 
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where the last term com es fro m the Ao = 0 eigenvalue. The stationary solution to the 
Fokker-Planck equation (4.1) is given by writing Q(cr) = e -V ( CT ). Doing so we reach 
the equation for the force 

dV{a 


d\ v 


= - 2 - 


M v 

'M vl 




4 A, 


, \l-\l 


Integrating over A„ and summing over v we obtain for the potential 

V (<t) = E (^ - 2 ln(A,)) - 2 Y, ^ I A 2 , - A^ |, 


which equates to the following JPDF for the eigenvalues. 

Q(a) = C N Y[\l n (j2 X l~ d »/ N 


n<m 


(44) 


This corresponds to the fixed trace symmetric GUE ensemble (see Section 3.4 of !)■ 


7. Real Wishart Bernoulli ensemble 


As a final example, we turn our attention to Bernoulli matrices B of size N x M 
(we shall assume N > M without loss of generality). Once again, every element is 
chosen independently from the set {±1 }/VN with equal probability and we denote this 
ensemble Bn as before. Such matrices are related to ensembles of random undirected 
graphs, as outlined in point ([ii| of the introduction. In the case N ^ M the matrix B 
does not have any eigenvalues, however there exists a factorisation of the form 

B = UAV t 


where A is an N x M matrix of singular values s v on the diagonal and the columns 
of the N x N matrix U and M x M matrix V give the left and right singular vectors 
respectively. The eigenvalues A„ of the matrix 

W = B t B. (45) 


are the square of the singular values, i.e. A„ = with (real) eigenvectors S\v) = \„\v) 
provided by the columns of V. Due to the nature of B , W commands a fixed trace 
Tr(W) = M for all B £ B N . 

As in the previous two examples, our Bernoulli ensemble has a Gaussian 
equivalent, known as the Wishart-Laguerre ensemble - corresponding to matrices W 
when the elements of B are normally distributed. In this case the JPDF for the 
eigenvalues can be determined exactly for any N, M, either by orthogonal polynomial 
techniques 43 44 , or via an appropriately adapted version of Dyson’s Brownian 
motion approach 45-47 . Analysis has also been performed for the eigenvalues of 


fixed-trace Wishart-Laguerre ensembles (see e.g. 35 and references therein). 

We will be interested in the limit N —> oo, with the ratio M/N — »c<l. In this 
case it is known that the mean spectral density for our Bernoulli ensemble tends to 
the following distribution 


P( A) = 


1 y/(a+ - A)(A - al) 
2irc A 


(46) 


where a± = (1± \/c) 2 and is supported in the interval [a_, a+]. This is the Marclienko- 
Pastur law [42] and implies our mean level spacing is of order 0(N~ 1 ). We are 
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interested in the fine scale behaviour of the eigenvalues of matrices (451. Recently 


this has been addressed in 19 , where the authors develop techniques used to treat 


Wigner matrices 10 , in order to analyse matrices B with elements of quite general 


probability distribution. 

Once again, we want to set up a random walk through Bn, analogous to that in 
Section [ 5 J The matrices B also represent the corners of a d jv-dimensional hypercube 
but this time with d n = NM, the number of independent parameters in B. A 
Hamming distance of 1 on the hypercube therefore corresponds to swapping just one 
element B pq , which means the single step perturbation SB = B' — B is the rank 1 
matrix 


6&W = -2B pq \p)(q\. 


(47) 


Similarly, a Af-step perturbation A B, given by (11), leads to the following 
perturbation in W 

W + AW = (B + AB) t {B + A B) = W+ [B T AB + A B T B + AB T AB}. (48) 

The task is to evaluate the moments of the spectral evolution due to a change in the 
matrix of AW. The details of which shall be omitted once again, since they follow 
very closely what has been outlined already in Section [5j Like in the previous two 
cases errors will inevitably be incurred because of the finite time effects of our random 
walk. These can be estimated thanks to delocalisation estimates of the eigenvectors 


which are applicable in covariance matrices 19 39 and results concerning spectral 


rigidity 19 for Wigner matrices. 


7.1. Evaluation of moments 

We begin with the expectation E(AW) of the change in the matrix due to a random 
walk of Af-steps. Using the expansion in terms of A B given by (48) and the analogous 
expression (|25| we find 


At 




JV=0 


Y[B T SB U 


+ (SB Ur ) T B] + Y {SB^fSB^ 


r,sGui 


. (49) 


In contrast to the previous two cases this contains quadratic terms in the matrix 
elements of SB, however this only presents an additional source of error that can be 
calculated. The leading term comes from a combination of the linear terms and the 
diagonal terms in the quadratic part, this is 

N,M 

E(AW) = Arj Y iB T 5B pq + ( 5B pq ) T B + (6B pq ) T {5B pq )\ . 
p,q= 1 


Inserting the rank 1 matrix (471, summing over p and q and applying this to the first 
term in the perturbation formula ( |23| ) we get 

E«i/|A W\v)) 


A p 


- = {vW m -W)\v)+0(At() =4(1—A,)+O(iV 2c - 2 ).(50) 


The leading correction in fact comes from the estimation of PAt(Af), as performed 
in ( [29| . The off-diagonal terms in the quadratic part of (49) provide a correction of 
order 0{Arf) = <D(N C ~ 2 ) which is slightly smaller. 
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If we proceed in a similar manner then we also obtain the second order 
contribution to the perturbation formula 

E(I(HAW|//)|-) + 0(jV 2c- 2)) X2(„\sw pq \p) 2 . 

An 

p,q 

Removing the absolute value is allowed since the inner product in real. In the previous 
two sections the analogous quantity obtained by summing over all matrix elements 
could be calculated exactly, whereas here we will will only get our desired expression 
up to some error. This is because SW itself is comprised of quadratic terms in SB, 
leading to cubic and quartic terms when SW is squared, unlike the previous cases. 
Writing this out explicitly we find 

5>l 6W pq \p.) 2 = \{ V \B T 8B pq \p) 2 + {v\{SB pq ) T B\p) 2 
p,q va 

+ 2{v\B T SB pq \n){v\{5B pq ) T B\p) 

+ (iy\(6B pq ) T 5B pq \p) 2 + 2(v\(SB pq ) T B\ri(v\(5B pq ) T 6B pq \p) 

+ 2(v\B T 5B pq \n)(p\(8B pq ) T SB pq \pL)]. (51) 

Inserting the expression ( |47| ) and summing over p and q we find 

Y / (v\tW pq \p,) 2 = ^[\„ + \ li + 25^\ u ) + 0(N- 2 ), (52) 

p,q 


where the leading order expression comes from the first three terms, quadratic in SB, 
and the remaining correction comes from the final three terms 

Yvy = o{N- 2 ). 


-) (M-A„-A m ) 


Here we have used that ^2 q n 2 p 2 « M/M 2 , 
happens with a high probability [l9j|39 . Collating (p50|) and ([52 1 


if they are delocalised, which we know 
together with the 


perturbation formula (23) we obtain the leading order expressions for the first two 
moments 




— 4 (1 — A„) + — Y 


+ Ay 


N Y Ay - A^ 


and 


My 


- 16 a 


(53) 


Again, the main corrections to these will come from higher order terms in the 
perturbation formula, which will rely on estimates regarding the rigidity of eigenvalues, 
recently obtained via local semicircle arguments 19 . One may also verify the off- 
diagonal, second order moments M ufi are of a sufficiently lower order than the diagonal 
ones Myy and can be neglected in the large N limit. Thus if the moments (53) are 


inserted into the Fokker-Planck equation the one may verify the resulting stationary 
solution is given by 


Q(<T)=n A 


i(JV-M+l)-l 


11 |Ay-A M |<5 ^Ay-M 


V<[1 


(using Tr(S') = M), which is the fixed trace JPDF for the Wishart/Laguerre ensemble. 
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In summary, we have shown how, in the limit of large matrix size, discrete random 
walks through the space of various Bernoulli ensembles lead to a stochastic motion that 
converges to the well-known Dyson Brownian motion model for the classical Gaussian 
ensembles. The key difference in our approach is that we do not assume any Brownian 
motion of the matrix elements from the outset, instead this motion arises simply from 
changing the sign of the matrix elements, randomly one by one. By calculating the 
moments of this evolution we obtain an approximate Fokker-Planck equation (19), 
and show the error decreases as the matrix dimension becomes larger. If one were to 
neglect this error then stationary solution would correspond to JPDF of certain well 
known fixed trace ensembles with Gaussian elements. 

We also believe that the methods could be adapted to treat ensembles in which 
the matrix elements are chosen from probability distributions other than Bernoulli. If 
these distributions are sufficiently well-behaved then again, we would expect similar 
results based upon the ideas of the CLT mentioned in the paper. 

During this work we have relied on the results of previous authors regarding 
the delocalisation of eigenvectors and the rigidity, or gaps, in eigenvalues. These 
were necessary in order to justify neglecting higher order terms in the perturbation 
theory, as well as higher moments stemming from the Taylor expansion ©• Since 
we would like to focus on the eigenvalues, it would be interesting to see whether 
such necessities can be circumnavigated. In fact, in a subsequent paper [49] we 
show how this can be achieved for similar unconstrained Bernoulli ensembles, as 
appear here. The approach is to develop a Fokker-Planck description for the traces 
tn = EllAU, 1 < n < N and then perform a transformation of variables back to the 
eigenvalues. The advantage of this method is that, in contrast to the approach in the 
present paper, it is non-perturbative. It may therefore offer a way to treat constrained 
random graph ensembles such as the random regular tournaments and random regular 
graphs outlined in the introduction, for which presently current methods cannot be 
applied. 

The methods may also be applicable to the spectral statistics of other graph 
operators, such as the graph Laplacian L = —A + D (where D is a diagonal matrix of 
the vertex degrees di = JA A;.?) or discrete Schrodinger operators H = L + V (where 
V is a random diagonal potential). These offer interesting alternatives to the matrices 
considered here, since one observes differences in the spectral properties, such as the 


mean spectral density [50 or the localisation of eigenvectors 51 52 
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Appendix A. Higher moments 

In the following we estimate the leading contributions coming from the third and 
fourth moments. We do this because these are odd and even moments and so their 
contributions will be slightly different. All others are of a similar suit and therefore 
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Appendix A.l. Cubic terms 


The linear and quadratic components provide us with the leading order contributions 
to the first and second moments, however we must also estimate higher moments. This 
means calculating the leading contributions to the moments M v)lK , etc., which 

will allow us to provide error estimates in our Fokker-Planck equation. In addition, 
as we shall discover, the main correction to M v will come, not from the errors arising 
in the calculations of the linear and quadratic terms, but from contributions coming 
from higher orders in the perturbation expansion (231. The third term of which is 
given by 


EE 


(is | AB | p) (p AB | k) (k \ AB \ is) 


(A,-A m )(A v -A k ) 'E (A.-A^) 2 

This requires us to estimate quantities of the form 


IHAb|m)I 2 HA.b|i/) 


(A.l) 


E ((v | AB | p) (p | AB | k) (k | AB \ is)). 


(A.2) 

Once again, we must follow the same procedure as before by summing over paths 
w, however the split into diagonal and off-diagonal terms is not quite so simple. The 
basic idea goes back to the fact that we have a sum over (almost independent) random 
variables, A B u = 5B Ur and that terms of the type E(|(i/|<5.B|^) | 2 ) will give a larger 
contribution than the \Ei((v\8B\p))\ 2 . Therefore when v = p for instance we have that 


(A.2) is approximately equal to 

At 2 E«i/|<SB|i/))E(|<i/|$B|K> 


| 2 ) + corrections, 

which gives a leading contribution ArfO(N c ~ 3 ) (if all three is, p and k were different 
then it would be 0). However we know from [lj the denominator in ( |A.l ) is of order 
0(N~ 2 ~ q ) with a probability that tends to 1 in the limit of large N (see Section [ 4 ] 
Therefore we get a correction to of order 0{N c+q ~ 1 ), which is in fact the dominant 
error to this time-averaged moment. A more detailed calculation of ( |A~2| begins with 
summing over all relevant paths and expanding A B u in terms of components uy £ w, 
just as we have done before in Section [5.1.1| and Section 5.1.2 


i.e. 


^ Pa* (A) 

x=o 


l^xl 


E E m 6 BuJt i«)(»i SB"*\v). 

wGfix r,s,u 

This will give us contributions coming from paths uj that contain either one, two or 
three elements that are the same. So we will get terms of the type 
At 


Pa* (A) 


x=o 1 1 

where 


<f> 


E/( 


e,e,e) + 


ci,e 2 


f(e i,ei,e 2 ) + $x /( e i 7 e 2,e 3 ) 


/(e i,e 2 ,e 3 ) = (is\SB ei \p)(p\8B e2 \K)(n\SB e3 \is) 

and as before e, £ { pq} P < q ■ We emphasise this is not exact but suffices to illustrate how 
the main contributions arise. In particular (comparing to equation (33)) the primed 
summands represent summations in which all the variables are not equal to each other. 
Though removing this restriction will only produce negligible errors. Moreover, we will 
have additional terms of the form /(e 1 , e 2 , ei), /(e 2 , e\, ei) etc. but these will be of the 
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same order as equivalent terms (i.e. /(ei, ei, ef)) and so since these are finite in number 
will only lead to some constant, independent of N, multiplying the contribution from 
/(ei,ei,e 2 ), which we neglect anyway. The three terms give respective contributions 

0(Ar]) Y Bp q {vpPq + is q p p ){p p n q + K q pbp)(npV q + n q v p ). (A.3) 

PA 


<W) £ Bpq(l/p[Aq H - Vqjlp) j ( ^ ^ Bpq([lpKq + ftqfJ>p)(Kpl / q KqlSp) 

\PA J \ pa / 

and 


(A.4) 


<W) Y, BpqijSpflq + ISqflp) J I ^ ^ Bpq([ApK,q Hq/lp ) J I ^ ^ Bpq^KplSq tiqVp) J .(A.5) 


, PA 


, pa 


< pa 


where again we have replaced the summation Y2 p < q with ^2 P q at the expense of a 
negligible error. Let us now estimate each term. For ( |A.3 ), if we have two equal 
eigenvalues, e.g. v — /x, then it becomes 

At jO{N~ l ) ^ 2 BpqVplSq(lSpK q + VqKp) 1 « Af]0(N~ 3 ) ^ B, 


pq\Bp^q T ft’pKq]’, 


p,q p,q 

where we have assumed again that the eigenvectors are delocalised, meaning that we 
take is 2 and k 2 = 0(N~ 1 ). Since the summation is of order unity, we can estimate 
(A.3) to be 0(N~ 3 )Arj. From the third term in the perturbation expansion (A.l) we 


see that at most two of is, p and n may be the same. If all three are different then 


(A.4) is 0 due to the first summation, however if we have is = p then (A.4) is equal to 


0(A V 2 )^ = 0{N c ~ 3 )Arj. 


Finally one can see the contribution ( |A.5 ) is 0 (up to subleading corrections of course). 
Thus we see that the dominant correction comes from (A.4), which is what one would 


expect from the central limit theorem, as it is quadratic in the time difference At. 


We now turn our attention to quantities of the form 
E ((is | AB | is) (p | AB | p) (k | AB \ k) ). 

These provide the leading order contributions to the third order moments 
and are of a slightly different to (A.2). Nevertheless they also give the same order 

abtain 

(A.6) 


contribution in N. To see this we proceed in exactly the same as before to obtain 
contributions of the type 


0{Arf) Y Bl q lSplS q PpP q KpK q 


o(A v 2 )IYb p 


£3 


pq[lp[iqKpKq 


, PA 


, PA 


and 


0( Av 3 ) Y BpqlSplSq £ BpqfJipflq £ BpqKpKq 


(A.7) 


(A.8) 
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Once again, the linear time contribution (A.6) is a fluctuating quantity that we 
estimate to be very small. The cubic contribution is given by A , qO{N 2c ~ 4 )X v X ll X Kj and 


is subdominant. This leaves (A.7). However, as we have already seen in Section 5.1.2 
this involves Y^ p q B 2 q p p pL q K, p n q = S^/N. Therefore if p = k this gives a contribution 
0(N C ~ 3 ) A?? and so applies for all terms of the type in which two or more of 

the eigenvalues are the same. 


Appendix A.2. Quartic terms 

The evaluation of quartic terms follows exactly the same line of reasoning. They 
appear in two forms, either 

E((^|AB| Ai )( / ,|AH| K )(K,|AH|7)( 7 |AB|^)), (A.9) 


which appear in the higher order terms of the perturbation expansion (23), or as the 
leading contributions to the moments M vpiK1 , for which they take the form 

e ((v|AH|i/)(/x|AB|/x)(k|AB|«){7|AH|7)). 

Let us briefly outline the latter case. Suppose that all parameters are equal, as is the 
case for the moments M vvvv . 

(A.8) then we see that any contributions with a factor „ B 


If we compare with the expressions (A.6), (A.7) and 

^ ^pqV p V(f) — Xy/dj^ 




will give a subleading contribution. In other words we have that ~Ei((v\8B\v) 2 ) is of a 
higher order in N than ~Ei((y\8 B\v)) 2 . This leaves us with two potential contributions. 
The first comes from the diagonal term, which is 


0(Ar,) = A V 0(N~ 4 ). 


P:Q 


The second comes from 


0{Ar, 2 )[Y j B : 


2 2 2 
pqVpV q 


= A pO{N c ~ 4 ). 


Thus we see the latter is of a larger order in N since c > 0. This will also be the case 
when we have two pairs of equal eigenvalues, i.e. moments of the type are also 

of the same order, however this is not the case when at least 3 eigenvalues are different 
and so moments of the type for instance are sub-leading. 

In addition, we may also find terms of the first type (A.9) which exhibit 
also contribution of order Ar]0(N c ~ 4 ). This means they provide an error to M„ 
approximately 0(N c+3s ~ 1 ), since terms in the denominator of the fourth term in the 
perturbation expansion are of the type (A„ — A M ) 3 = 0(N~ 3 ~ 3q ). Therefore we see 
that fourth order terms give a marginally larger contribution to M v than third order 
terms, which contribute 0(N c+2q ~ 1 ). However, beyond this we find 0(N 2c+4q ~ 2 ) and 
0(N 2c+5q ~ 2 ) for the fifth and sixth orders, which are subleading in comparison to the 
fourth. 
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